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1 . 0  INTRODUCTION 


This  document  is  Volume  II  of  a  two-volume  Final  Technical 
Report  which  summarizes  a  major  task  carried  out  by  the  Scientific 
Studies  Corporation  (SSC)  in  Phase  II  of  the  "Multichannel  System 
Identification  and  Detection  Using  Output  Data  Techniques"  Small 
Business  Innovation  Research  (SBIR)  program  for  Rome  Laboratory 
(RL) .  The  companion  Volume  I  presents  the  formulations,  analyses, 
and  simulation  results  obtained  in  the  program  (Roman  and  Davis, 
1996).  Distribution  of  Volume  I  is  limited  as  dictated  by  the 
standard  SBIR  regulations,  whereas  distribution  of  this  volume  is 
unlimited. 

This  volume  discusses  the  analytic  and  software  (MATLAB- 
based)  model  for  the  multichannel  output  waveform  in  airborne 
surveillance  phased  array  radar  systems.  Specifically,  target, 
jamming  (spatially-localized,  broadband  interference),  receiver 
noise,  and  ground  clutter  components  are  modeled.  Of  these 
components,  ground  clutter  is  the  most  difficult  to  model  in  the 
context  of  surveillance  array  radar.  The  ground  clutter  model 
follows  the  approach  used  by  Jaffer  et  al .  (1991),  and  more 
recently  by  Ward  (1994)  . 

A  description  of  the  SSC  model-based  multichannel  detection 
methodology  in  the  context  of  space/ time  processing  for  phased 
array  radar  systems  is  provided  in  Volume  I  (Roman  and  Davis, 
1996).  In  order  to  validate  the  methodology  and  to  evaluate 
alternative  model  identification  algorithms  in  the  context  of 
space/time  radar  applications  it  is  necessary  to  have  either  real 
or  simulated  radar  data.  Since  real  data  was  unavailable  to  this 
program,  SSC  developed  analytic  and  software  models  to  generate 
simulated  data.  For  the  purposes  of  this  program  it  suffices  that 
the  simulated  data  be  "representative"  of  real  data,  while  placing 
modest  requirements  on  computer  resources  (storage;  execution 
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time) .  The  software  model  should  be  flexible,  in  order  to  admit 
future  enhancements  in  a  straightforward  manner.  Additionally, 
the  model  must  include  a  representation  for  the  "true"  channel 
output  covariance  matrix  sequence  to  allow  determination  of 
performance  in  the  "known  covariance"  case.  This  latter  feature 
is  important  for  algorithm  evaluation,  and  is  unavailable  with 
real  data. 

The  conditions  and  assumptions  for  the  model  formulation  are 
(^0gQ2ribed  in  Section  2.0,  and  the  true  covariance  model  and  the 
data  generation  model  are  presented  for  each  component  (signal; 
jamming;  noise;  clutter) .  Simulation  results  are  presented  in 
Section  3.0  to  allow  comparison  with  corresponding  results 
presented  by  Ward  (1994) .  A  summary  and  suggestions  for  further 
work  are  presented  in  Section  4.0. 
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2 . 0  AIRBORNE  SURVEILLANCE  PHASED  ARRAY  RADAR  MODEL 

Airborne  surveillance  for  moving  target  detection  using 
phased  array  radar  is  one  of  the  major  thrusts  at  RL. 
Consequently,  the  scenarios  and  systems  considered  in  the  model 
defined  herein  reflect  that  significance.  The  model  is  defined 
for  a  side-looking  phased  array  radar  configuration,  with  scanning 
capability  of  ±90  deg  from  boresight.  Each  component  (target; 
noise;  jammer;  clutter)  is  modeled  independently  of  the  others. 
The  structure  selected  to  model  each  component  allows  generation 
of  the  true  covariance  matrix  sequence,  as  well  as  generation  of 
independent  statistical  data  realizations . 

Moving  targets  are  modeled  as  point  targets  in  azimuth  and 
elevation,  and  each  moving  target  is  modeled  as  the  output  of  a 
first-order  linear  system  in  state-space  form.  This  allows  for  a 
wide  range  of  conditions  (including  Swerling  Cases  1  through  4  as 
well  as  the  deterministic  case),  in  a  simple  and  general  format. 
With  this  structure  Gs  targets  are  grouped  as  a  Gsth-order  linear 
system.  Barrage-type  jammers  are  considered,  and  modeled  as  point 
sources  of  Gaussian-distributed  broadband  noise  in  azimuth  and 
elevation.  Receiver  noise  is  modeled  as  uncorrelated  in  time 
(pulse-to-pulse)  as  well  as  in  space  (channel-to-channel) .  Ground 
clutter  is  modeled  based  on  the  principle  of  a  large  number  of 
statistically-independent  clutter  patches,  and  each  clutter  patch 
consists  of  a  large  number  of  individual  re-radiators. 

2 . 1  Scenario  Description  and  Problem  Formulation 


Consider  a  linear  array  radar  consisting  of  J  equally-spaced, 
identical  antenna  elements  (or  identical  beamformed  sub-arrays)  in 
a  side-looking  configuration  on  a  surveillance  platform  moving  at 
a  constant  velocity  Vp  in  level  flight  at  an  altitude  Hp.  Such  a 
scenario  is  depicted  in  Figure  2-1.  As  suggested  in  Figure  2-1, 
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the  array  is  aligned  with  the  aircraft's  longitudinal  axis,  and 
the  aircraft  velocity  makes  a .positive  crab  angle  y  with  the 

aircraft's  longitudinal  axis.  For  transmission,  a  power  pattern 
G,((|),6)  is  assumed  for  the  array,  where  (])  is  the  azimuth  angle,  and 
0  is  the  elevation  angle.  The  array  main  beam  can  be  steered  in 
azimuth  and  elevation,  with  a  positive  main  beam  azimuth  angle  (t)o 
measured  clockwise  from  the  array  boresight  (a  positive-valued  (|)o 

is  shown  in  Figure  2-1) ,  and  a  positive  main  beam  elevation  angle 
6q  defined  as  a  rotation  up  from  the  local  horizontal  plane  at  the 

radar  array  location  (a  negative-valued  ©g  is  shown  in  Figure  2-1)  . 
For  simplicity,  the  elevation  dimension  is  ignored  herein,  and  the 
clutter  return  is  calculated  for  the  mainbeam  elevation  only. 

Thus,  the  antenna  power  pattern  used  to  generate  the  ground 
clutter  model  is  Gt((l),6o)-  The  radar  array  is  radiating  a  coherent 
pulse  train  at  a  constant  radiation  frequency  Iq  (with 
corresponding  radiation  wavelength  ^q)  ,  and  a  constant  pulse 
repetition  frequency  (PRF)  fpRp.  For  reception,  the  power  pattern 

of  each  individual  linear  array  element  (or  sub-array  beamformed 
into  a  linear  array)  is  denoted  as  gr((t)i0);  as  in  transmission,  the 

receive  power  pattern  is  restricted  to  the  mainbeam  elevation 
only,  0  =  0Q. 

It  is  assumed  that  the  range  to  moving  targets,  ground,  and 
to  broadband  interference  sources  is  large,  so  that  each  wavefront 
arriving  at  the  array  is  flat  (in  the  scale  of  X.^)  .  This  is  true 

for  most  surveillance  scenarios.  It  is  assumed  also  that  each  of 
the  received  wavefronts  is  statistically  stationary  (in  the  wide 
sense)  for  N  pulses  from  a  given  range  gate,  the  duration  of  one 
coherent  processing  interval  (CPI) .  The  4/3  Earth  approximation 

(Nathanson,  1991)  is  adopted  for  determination  of  the  grazing 
angle,  \|/(,.  Table  2-1  lists  all  variables  (physical;  scenario; 

system)  utilized  in  the  model;  this  table  also  includes  the  key 
parameters  associated  with  the  return  from  a  clutter  patch 
(normalized  spatial  and  Doppler  frequencies;  geometry  parameters) . 
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Figure  2-1.  Airborne  surveillance  scenario  and  geometry. 
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fs  =—  cos(e)  sin(<t))  =  0.5  cos(0)  sin((t))  d  =■ 


Normalized  Doppler  frequency: 
cos(0)  sin((|)- y)  v= 

fppp^C 

f  2Vr, 

fj  =  — = - E_  cos(0)  sin((|)  -y)  =  O.5-ucos(0)  sin((t)-Y) 

^PFF  fpRF^C 


Clutter  patch  parameters 


COS(Vc) 

.  J  r|  +  HD(HD  +  2rc) 

0„  =  -sin-i 

°  L  2rc(Hp  +  r£) 


Ad)  =— 22L 
N, 


V-  =  sin 


mainbeam  clutter  return  azimuth  angle  (shown  positive  in  Figure  2-1) 
target  return  azimuth  angle  (shown  positive  in  Figure  2-1) 
interference  source  azimuth  angle  (shown  negative  in  Figure  2-1) 
mainbeam  clutter  return  elevation  angle  (shown  negative  in  Figure  2-1) 
target  return  elevation  angle  (not  shown  in  Figure  2-1) 
interference  source  elevation  angle  (not  shown  in  Figure  2-1) 
aircraft  platform  crab  angle  (shown  positive  in  Figure  2-1) 
platform  velocity 
platform  altitude 

target  radial  velocity  (not  shown  in  Figure  2-1) 

array  inter-element  spacing 

radiation  wavelength 

Doppler  shift  frequency 

pulse  repetition  frequency 

clutter  ring  grazing  angle  (shown  positive  in  Figure  2-1) 

range  to  clutter  ring 

clutter  patch  range  extent 

clutter  patch  angular  extent 

radar  range  resolution 

number  of  clutter  patches  in  the  clutter  ring 

effective  mean  radius  of  the  Earth 
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In  most  surveillance  array  systems  the  narrow  bandwidth 
criterion  is  satisfied,  and  it  is  assumed  to  be  true  for  the  cases 
considered  herein.  For  a  J-element  linear  array  with  inter-element 
spacing  d,  the  narrow  bandwidth  criterion  can  be  expressed  as 
(Hudson,  1981) 

V 

(2-1)  (J-1)d« -  Narrow  Bandwidth  Criterion 

^BBW 


where  V  (lower  case  greek  nu)  denotes  the  speed  of  propagation  of 
the  electromagnetic  wave  and  f0BW  one-sided  amplitude 

waveform  (baseband)  bandwidth.  In  this  work  the  inter-element 
spacing  d  is  selected  as  (see  Table  2-1) 


(2-2) 


This  choice  of  inter-element  spacing  avoids  the  generation  of  more 
than  one  array  antenna  beam,  the  so-called  grating  lobes  (Monzingo 

and  Miller,  1980;  Skolnik,  1970).  Many  practical  radar  array 
systems  are  designed  with  d<0.5^Q,  which  also  results  in  one  beam. 

An  equivalent  expression  for  the  narrow  bandwidth  criterion  is 
obtained  by  substituting  Equation  (2-2)  into  Equation  (2-1), 

2fr 

(2-3)  J-1« - ^  Narrow  Bandwidth  Criterion 

^BBW 


This  criterion  may  not  be  met  for  arrays  with  a  large  number  of 
elements  and/or  a  large  amplitude  waveform  bandwidth. 

The  elevation  dimension  is  included  to  first-order  in  this 
formulation.  This  approximation  is  reasonable  for  long-range 
surveillance  radar,  and  suffices  for  the  purposes  of  this  program. 
A  more  complete  model  includes  a  planar  (two-dimensional)  array 
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with  main  beam  elevation  steering.  However,  the  model  adopted 
herein  can  be  used  to  represent  the  cases  where  the  channels 
correspond  to  antenna  subarrays  in  either  one  or  both  dimensions 
(azimuth  and  elevation);  that  is,  the  channels  are  J  subsets  of 
appropriately-combined  individual  elements.  In  particular,  the 
clutter  component  model  includes  the  cases  wherein  Je  elevation 
channels  are  beamformed  into  a  single  row. 

Consider  the  output  of  the  ith  channel  corresponding  to  a 
single  range  resolution  cell  after  pulse  compression,  demodulation 
down  to  baseband,  and  temporal  sampling.  Notice  that  for  a  single 
range  resolution  cell,  each  channel  output  includes  the  return 
from  multiple  ground  clutter  rings  (the  primary  range  ring  as  well 
as  ambiguous  range  rings) .  Denote  this  complex-valued,  discrete¬ 
time  waveform  as  {Xj(n)  I  n  =  0,  1,  .  .  .  ,  N-1}  =  {Xj(n)}.  Specifically,  for  a 
single  range  resolution  cell  and  one  CPI,  the  ith  channel  output  is 

(2-4)  Xj(n)  =  Sj(n)  +  ii(n)  +  Cj(n)  +  Wj(n)  0<n<N-1 

where  n  is  an  integer  denoting  the  sampling  time  instant,  {Sj(n)} 
denotes  the  signal  (moving  target)  return  at  the  ith  channel,  {ij(n)} 

denotes  the  total  broadband  interference  (jamming)  received  at  the 
ith  channel,  {Cj(n)}  denotes  the  total  ground  clutter  return  at  the 
ith  channel,  and  {Wj(n)}  denotes  the  combination  of  all  the 
statistically-independent  noise  sources  in  the  ith  channel.  Array 
receiver  noise  is  usually  the  dominant  noise  contributor  to  {Wj(n)}, 
and  that  is  assumed  to  be  true  for  the  model  considered  herein. 
Each  of  the  output  components  is  assumed  to  be  a  stationary, 
ergodic,  zero-mean,  Gaussian-distributed  random  process. 

The  temporal  sampling  interval  At  for  the  channel  output 
waveform  is  the  pulse  repetition  interval  (PRI)  T|,  which  is  the 

inverse  of  the  PRF.  Thus,  it  is  more  appropriate  to  express  the 
function  argument  in  Equation  (2-4)  as  nT|,  instead  of  n.  However, 
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a  PRF-normalized  frequency  variable  is  used  in  the  formulation 
developed  below,  and  such  normalization  is  equivalent  to 
normalization  of  the  sampling  interval  by  the  PRI .  As  a  result, 
the  sampling  interval  At  associated  with  the  variable  n  (or  with 

the  variable  m,  which  is  used  herein  to  denote  correlation  lags), 
can  be  assumed  to  be  equal  to  unity  and  dimensionless.  Also,  the 
initial  time  of  the  CPI  is  taken  as  n  =  0  without  loss  of  generality 
since  the  process  is  stationary. 

Let  channel  1  be  the  temporal  and  spatial  reference  for  the 
array.  For  a  uniform  linear  array  the  waveforms  from  the  J 
channels  can  be  concatenated  into  the  following  vectors, 

(2-5)  .  x(n)  =  s(n)  +  i(n)  +  c(n)  +  w(n)  0<n<N-1 

"xi(n) 

X2(n) 

(2-6)  x(n)  = 

-Xj{n) 

"si(n)1  r  Si(n) 

S2(n)  s^(n+T) 

(2-7)  s{n)  =  .  = 

.Sj(n)J  L  s.|(n+Ki:) 

'ii(n)1  r  ii(n) 

i2(n)  ii(n+T) 

(2-8)  i(n)  =  .  = 

Jj(n)J  L.ii(n+Kx) 
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Ci(n) 

C2(n) 

(2-9)  c(n)  = 

-Cj(n) 

"  Wi(n) 

W2(n) 

(2-10)  w(n)  = 

.  Wj(n) 

with  K  =  J-1  in  order  to  simplify  notation,  and  where  t  is  the 

inter-element  time  advance  (with  respect  to  channel  1,  the 
reference  channel)  .  For  a  range  resolution  cell  at  azimuth  angle  (j) 
and  elevation  angle  0,  the  inter-element  time  advance  is  obtained 
from  the  geometry  established  in  Figure  2-1  as 

d  d 

(2-11)  X  =  — cos(0)sin((|))  = - cos(0)sin((l)) 

V  Xcfc 

where  all  variables  are  as  de,fined  previously  (see  Table  2-1)  . 

As  indicated  in  Equations  (2-7) -(2-9),  the  kth  element  of 

each  of  the  signal,  interference,  and  clutter  components,  is  an 
advanced  version  of  the  scalar  waveform  in  channel  1,  with  T  as  the 

inter-element  advance.  This  follows  from  the  assumptions  and 
conventions  established  above.  Array  noise,  however,  is  assumed 
to  be  uncorrelated  from  channel  to  channel,  and  with  the  same 
power  in  each  channel  (the  equal  power  assumption  can  be  relaxed 
easily,  but  it  is  representative  of  realistic  conditions) .  These 
characteristics  determine  the  structure  of  the  covariance  matrix 
sequence  of  each  component. 


Ci(n) 

(n+x) 

(n+Kx) 
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2.2  Signal  (Moving  Target)  True  Covariance  Model 


The  signal  (moving  target)  component  at  the  output  of  the  kth 
channel,  {Sk(n)},  is  the  radar  return  from  each  moving  target  located 

within  the  range  resolution  cell  of  interest.  Consider  a  target 
moving  towards  the  radar  array  with  a  constant  velocity  Vj  along  a 
trajectory  aligned  with  azimuth  angle  (j)j  and  elevation  angle  6j. 

The  target's  radial  velocity  is  defined  to  be  positive  when  the 
target  is  moving  towards  the  radar,  and  negative  otherwise.  Such 
a  target  has  a  normalized  Doppler  frequency  shift  fj^j  of  the  form 

(Figure  2-1  and  Table  2-1) 


(2-12) 


tD 


2[V^+Vp  cos(e^)sin((|)^-Y)] 


PRF 


^PRF^C 


where  is  the  normalized  target  Doppler  frequency  shift,  and  the 

PRF  provides  the  normalization.  This  normalized  Doppler  frequency 
assumes  unambiguous  values  in  the  range  -0.5  <  <  0.5 . 


From  the  definition  of  spatial  frequency  in  Table  2-1  and  the 

geometry  depicted  in  Figure  2-1,  the  target's  dimensionless 
spatial  frequency  fjg  is  of  the  form 

(2-13)  f,3  =fct  =  — cos(0,)sin((t)j) 

Kq 


where  the  inter-element  time  advance  x  is  as  defined  in  Equation 
(2-11) .  Adopting  the  half-wavelength  value  for  d  (see  Equation 
(2-2)),  the  target's  spatial  frequency  becomes 

(2-14)  =  0.5  cos(et)  sin((t),)  -90°  <  0,  <  90° ;  -90°  <  (|),  <  90° 

Given  the  allowable  range  of  values  for  the  azimuth  and  elevation 
angles,  it  follows  that  -0.5  <  fjg  <  0.5 .  Notice  that  the  sign  of  the 


11 


azimuth  angle  determines  the  sign  of  the  spatial  frequency.  For 
the  geometry,  conventions,  and  definitions  established  in  Figure 
2-1,  Table  2-1,  and  above,  the  return  from  a  target  illuminated  by 
the  main  beam  at  a  positive  azimuth  angle  has  a  positive  spatial 
frequency  fjg. 

Moving  targets  are  modeled  as  the  output  sequence  of  a 
linear,  stationary,  discrete- time ,  complex-valued,  stochastic 
system  in  state-space  form.  Specifically,  for  the  case  of  Gg 

targets,  the  model  adopted  herein  is 

(2-15a)  y(n  +  1)  =  Fy(n)  +  u(n)  n>0 

(2-15b)  s(n)  =  H^y(n)  n>0 

{2-15C)  E[y(0)]  =  0 

{2-15d)  E[y(0)yH(0)]  =  P(0)  =  P 

(2-15e)  E[y(n)yH(n)]  =  P(n) 

where  y(n)  is  the  Gg-dimensional  state  vector,  with  y(0)  Gaussian- 
distributed  and  zero-mean;  u(n)  is  the  Gg-dimensional  input  vector; 
and  s(n)  is  the  Gs-dimensional  output  vector.  The  input  sequence 
{u(n)}  is  Gaussian-distributed  with  mean  zero  and  covariance  matrix 

(2-16)  E[u(n)u^(n)]  =  Q 

It  follows  that  the  state  and  output  sequences  are  both  Gaussian- 
distributed  with  mean  zero  also.  Matrix  F  is  the  GgXGg  system 
matrix,  which  defines  the  time  evolution  of  the  process.  Matrix  H 
is  the  GgXJ  output  distribution  matrix,  P(n)  is  the  state 
covariance  matrix  at  time  n,  and  P  is  the  state  covariance  matrix 
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at  steady-state.  For  a  system  of  the  type  (2-15) ,  stationarity 
implies  asymptotic  stability,  which  in  turn  implies  that  steady- 
state  conditions  are  attained  at  a  finite  time.  The  system 
matrix,  output  matrix,  input  covariance  matrix,  and  state 
covariance  matrix  are  defined  as  follows: 


(2-17a)  F  = 


F(1.1)  0 

0  F(2,2) 


0 

0 

F(Gs,Gs) 


(2 -17b)  F(k,k)  =  Fk  =  exp(j2jiftdk) 


(2-17c)  0<ak<1 


(2-18a)  H  §j(^ts2)  ■"  -j(^tsGs)] 


(2-18b)  §j(^ts)  = 


1 

j2Kf,sK 


2-18C)  K  =  J-1 


(2-19a)  Q  = 


Q{1.1)  0 

0  Q(2,2) 

0  0 


0 

0 

Q(Gs,Gg) 


(2-i9b)  Q(k,k)  =  Qk  =  (1  -  a^)  Pr  =  (1 ' 4)  4k 


k  =  1 , 2, . . . ,  Gg 
k=  1, 2 . Gg 


k=1,2,  ...,G, 


(2-20)  P(n)  =i-  P 


as  n  oo 
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- 1 

_JD 

O 

O 

1 

0  •••  0 

(2-21)  P  = 

0  P2  •••  0 

= 

0  a^2  •••  0 

0  0  •••  Pp 

L  Gs  J 

0  0  -  <y^G^_ 

Selection  of  the  initial  state  covariance  matrix  as  a  diagonal 
matrix  equal  to  the  steady-state  covariance  (see  Equation  (2-15d) ) 
results  in  the  straightforward  calculation  of  the  the  input  noise 
covariance  (Equation  (2-19b) ) ,  and  simplifies  data  generation 
(Section  2.7).  Furthermore,  each  target  can  be  assigned  a 
different  power  since  the  diagonal  elements  of  P  can  be  distinct 
from  each  other.  This  allows  assignation  of  a  distinct  signal-to- 
noise  ratio  (SNR)  to  each  target. 

This  model  has  several  features  of  interest.  First,  the 
temporal  evolution  of  the  kth  target  is  represented  by  the  kth 
state  exclusively  via  the  damping  coefficient,  a,^,  and  the  target's 
normalized  Doppler  shift  frequency,  Second,  the  system  is 
exponentially  stable  because  all  the  coefficients  3)^  are  less  than 

unity,  as  defined  in  Equation  (2-17c)  .  As  a  consequence,  the 
stochastic  identification  algorithms  of  interest  in  this  program 
can  be  applied  without  concern  for  pathological  behaviour.  Third, 
the  spatial  structure  is  assigned  via  the  columns  of  H*^, 
separately  from  the  temporal  structure.  The  columns  of  defined 
in  Equation  (2-18b)  are  referred  to  as  spatial  frequency  vectors, 
where  fjg|^  is  the  spatial  frequency  of  the  kth  target.  Fourth,  the 

power  of  the  kth  moving  target  is  P,^  =  CJgi^ ,  and  this  power  is 

contributed  to  each  channel.  As  a  result,  the  total  signal  power 
in  each  channel  is  equal  to  the  sum  of  the  diagonal  of  P.  Fifth, 
the  covariance  matrix  sequence  of  the  output  is  determined  in 
closed  form  as  a  function  of  the  system  parameters  defined  above. 
Specifically,  let  {RgsCiTi)}  denote  the  JxJ  auto-covariance  matrix 
sequence  for  the  zero-mean  process  {s(n)},  then 
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(2-22) 


R3s(m)  =  E[s(n)sH(n-m)]  =  H'^F'^PH 


m  >0 


with  negative  lags  defined  as 

(2-23 )  RssC-iT^)  =  E[s(n)s'^(n+m)]  =  Rss(m)  m  >  0 

Equations  (2-22)  and  (2-23)  define  the  true  covariance  matrix 
sequence,  {RasCf^)}-  Matrix  Rss(fh)  is  Hermitian  only  for  m  =  0; 
however,  matrix  Rss(fh)  Toeplitz  structure  for  all  m  due  to  the 

spatially-induced  relationship  among  the  elements  of  s(n),  as 
expressed  in  Equation  (2-7)  . 

Equation  (2-22)  can  be  used  to  generate  the  true  covariance 
matrix  sequence  directly.  However,  due  to  the  special  structure 
of  matrices  F,  H,  and  P,  the  right-hand- side  of  Equation  (2-22) 
can  be  expressed  in  a  form  that  requires  less  computations. 
Namely,  the  mth  lag  covariance  matrix  can  be  expressed  as 

(2-24)  R.sCn)  =  i  a™  P,  e;(U  m  >  0 

k=1 


Equation  (2-24)  requires  less  computations  than  Equation  (2-22)  . 
However,  given  the  structure  of  the  spatial  frequency  vectors,  the 
form  of  Equation  (2-24)  ,  and  the  fact  that  matrix  Rss(^) 

Toeplitz,  a  simpler  approach  can  be  defined  to  generate  the  true 
matrix  covariance  sequence,  as  described  next.  First  compute  the 
following  two  J-dimensional  vectors, 

Gs 

(2-25)  =  X  ^k  (^tsk)  rn  >  0 

k=1 
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(2-26)  eH(m)  =  X 


m>0 


k=1 


Notice  that  ej,(m)  is  a  coliamn  vector,  whereas  ej^(m)  is  a  row  vector. 
Then,  for  each  desired  value  of  m,  matrix  Rss(fl')  is  formed  as 
follows:  (a)  set  the  first  column  of  Rss(fn)  equal  to  vector  §^,(01), 

which  defines  the  main  diagonal  and  all  the  diagonals  below  the 
main  diagonal;  and  (b)  set  the  first  row  of  Rss(m)  equal  to  vector 
e»,  which  defines  the  main  diagonal  and  all  the  diagonals  above 

the  main  diagonal.  There  is  no  conflict  along  the  main  diagonal 
because  the  first  element  of  ^(m)  is  the  same  as  the  first  element 

of  e[^(m). 

2 . 3  Broadband  Interference  (Jamming)  True  Covariance  Model 

The  broadband  interference  (jamming)  at  the  output  of  the  kth 
channel,  {ik(ri)},  represents  the  received  interference  from  all 
direct  jamming  sources.  Jammer  multipath  and  distributed  jammers 
are  not  included  in  the  model  at  the  present  time.  Consider  a 
jamming  point  source  located  at  azimuth  angle  (|)j,  as  shown  in 
Figure  2-1  for  a  negative-valued  angle  case,  and  elevation  angle  0j 
(not  shown  in  Figure  2-1) .  This  interference  is  modeled  at  the 
output  of  the  each  channel  as  a  Gaussian-distributed,  zero-mean, 
white  sequence  with  variance  cf .  The  spatial  variation  of  this 
component  is  described  by  the  dimensionless  spatial  frequency  fjs, 
which  is  given  as 

d 

(2-27)  f|s  =  fct  =  -7 — cos(ej)sin((|)j)  =  0.5  cos(0j)sin((t).) 

A-c 

where  the  inter-element  time  advance  x  is  as  defined  in  Equation 
(2-11) ,  and  the  half -wavelength  inter-element  spacing  assumption 
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(Equation  (2-2))  has  been  applied.  Each  channel  receives  the  same 
amount  of  power  from  this  jammer. 

The  true  covariance  matrix  sequence  {Rjj(m)}  for  G|  interference 
sources  as  described  above  is  given  as 

(2-28a)  R„(0)  =  Eli(n)iH(n-m)]  =  f 

k=1 

(2 -2  8b)  Rii(m)  =  Ea(n)iH(n-m)]  =  [0]  m  0 

where  ej(fjsk)  is  the  spatial  frequency  vector  for  the  kth  jamming 
source  at  the  spatial  frequency  fjgk,  and  is  the  variance  (power) 
of  the  kth  jamming  source.  Each  channel  receives  the  same  amount 

of  power  from  this  jammer.  Notice  that  Equation  (2-28a)  is  of  the 
same  form  as  Equation  (2-24)  for  m  =  0.  Thus,  matrix  Rjj(O)  is 
Toeplitz.  Also,  the  generation  of  Rjj(O)  can  be  carried  out  via  the 
same  procedure  used  to  generate  Rss(O)  • 

2  •  4  Receiver  Array  Noise  True  Covariance  Model 


The  receiver  array  noise  at  the  output  of  the  kth  channel, 
{Wk(n)},  represents  all  noise  sources  that  arise  within  the  kth 

channel.  This  component  is  modeled  as  a  Gaussian-distributed, 
zero-mean,  white  sequence  with  variance  (power)  a^k  /  and 

uncorrelated  from  channel  to  channel.  Thus,  the  spatial  variation 
is  broadband  also.  It  follows  that  the  true  covariance  matrix 
sequence  {Rww(lT^)}  this  output  component  is  given  as 


(2-29a) 


Rww(0)  =  E[w(n)wH(n)j^ 


0 

0 
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(2-29b)  Rww(m)  =  E[w(n)  w'^(n-m)]  =  [0] 


m  0 


Notice  that  the  covariance  matrix  Ryvyy(O)  as  defined  in  Equation  (2 
29a)  is  not  Toeplitz.  This  general  model  is  applicable  in  some 
cases,  specially  when  each  channel  corresponds  to  a  different  type 
of  sensor.  However,  in  the  surveillance  radar  case  considered 
xs  appropriate  to  assign  the  same  variance  to  the  noise  in 
each  channel.  Thus,  the  model  adopted  herein  for  the  receiver 
gxray  noise  true  covariance  matrix  sequence  {Rvwv(^)} 

{2-30a)  Rww(O)  =  E[w(n)  w^(n)]  =  Ij 

(2-30b)  Rww(m)  =  E[w(n)  w^(n-m)]  =  [0] 

^xtih  the  noise  sequence  variance  (power)  in  each  channel. 

Matrix  Rww(O)  in  Equation  (2-30a)  is  Toeplitz.  Notice  that  the 
difference  between  the  noise  covariance  Rww(O)  ^nd  the  interference 
covariance  R|j(0)  is  the  spatial  correlation  present  in  Rjj(O).  This 
difference  is  sufficient  for  adaptive  cancelation  algorithms  to 
distinguish  between  noise  and  interference. 

2 . 5  Ground  Clutter  True  Covariance — Model 

The  ground  clutter  component  at  the  output  of  the  kth 
channel,  {c,^(n)},  is  the  combination  of  the  radar  return  from  all 

ground  patches  irradiated  by  the  transmitted  pulse  sequence.  This 
includes  the  return  from  ground  patches  in  the  primary  range  ring 
(see  Figure  2-1) ,  which  is  irradiated  by  the  array  main  beam  as 
well  as  by  the  sidelobes;  also  included  is  the  return  from  range 
rings  at  ambiguous  ranges  and  from  ranges  illuminated  by  the 
elevation  pattern  sidelobes .  Since  each  range  ring  extends  over 
all  azimuth  angles,  ground  clutter  return  exhibits  power  in  a 
continuum  of  Doppler  and  spatial  frequencies.  In  fact,  a  typical 
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ground  clutter  power  spectrum  (see,  for  example,  Jaffer  et  al . , 
1991)  exhibits  a  "ridge"  in  the  two-dimensional  (2-D)  frequency- 
plane,  as  indicated  in  Figure  2-2.  In  cases  where  the  crab  angle 
is  zero,  the  clutter  ridge  occurs  along  a  straight  line  that 
passes  through  the  origin  of  the  2-D  frequency  plane,  as  in  Figure 
2-2.  When  the  crab  angle  is  non-zero,  the  trace  defined  by  the 
ridge  is  a  nonlinear  function  of  various  system  and  scenario 
parameters.  The  number  of  clutter  ridges  is  doubled  when  array 
backlobes  are  significant  and  the  crab  angle  is  non-zero  (Ward, 
1994) .  The  model  defined  herein  accounts  for  sidelobe  clutter  and 
backlobe  clutter,  but  does  not  include  ambiguous  range  rings  nor 
does  it  include  clutter  from  elevation  pattern  sidelobes. 

^cc(^d  ’^s  ) 


Figure  2-2.  Two-dimensional  clutter  spectrum  and  angle  a. 


The  clutter  ridge  angle  a  (see  Figure  2-2)  is  defined  as 
(Jaffer  et  al.,  1991;  Ward,  1994) 
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where  X>  (the  Greek  letter  upsilon)  is  the  slope  of  the  clutter 
ridge  (Table  2-1).  When  u >  1  the  clutter  spectrum  has  ambiguities 
(components  at  normalized  Doppler  frequencies  greater  than  0.5), 

and,  consequently,  the  clutter  Doppler  spectrum  aliases  into  the 
unambiguous  Doppler  space,  -0.5  <  <  0.5  (Ward,  1994)  .  This 

condition  is  manifested  in  the  2-D  frequency-domain  plane  as 
multiple  ridges  (straight  lines)  of  different  length.  If  the  crab 
angle  is  non-zero,  the  straight-lined  ridges  become  curved. 

The  ground  clutter  model  described  herein  is  based  on  the 
analytic  model  proposed  by  Jaffer  et  al .  (1991),  and  also  adopted 
recently  by  Ward  (1994).  However,  the  notation  used  herein 
follows  Ward  (1994)  more  closely.  Two  non-trivial  simplifications 
to  these  models  are  introduced  herein:  (a)  only  a  single  range 
ring  is  modeled,  and  (b)  intrinsic  clutter  motion  is  not  modeled. 
In  Ward's  formulation  the  clutter  return  from  a  range  ring  is 
determined  by  equal-area  patches,  and  each  patch  consists  of  a 

single,  complex-valued  scatterer  with  mean  zero  and  Rayleigh- 
distributed  magnitude.  Furthermore,  the  scatterers  are  assumed  to 
be  pairwise  statistically  independent.  Ward's  model  for  the  true 
covariance  sequence  and  for  data  generation  are  based  on  an 
expression  for  clutter-to-noise  ratio  (CNR)  as  a  function  of 
azimuth,  at  a  fixed  elevation  angle.  Specifically,  for  a  single 
clutter  ring  Ward's  CNR  as  a  function  of  azimuth  and  elevation  is 
of  the  form 
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(2-33) 


y<^.e)  = 


P|Tu4CoG„S,('t'.e)g,(4'.e) 

(4)tfN„L,r‘ 


where  Pj  denotes  the  peak  transmitted  power,  T^,  denotes  the 
uncompressed  pulse  duration,  Xq  denotes  the  radiation  wavelength, 
denotes  the  clutter  patch  radar  cross  section  (RCS) ,  denotes 
the  receive  column  pattern  {beamforming)  gain,  Gt((l),6)  denotes  the 

transmit  antenna  power  pattern  as  a  function  of  azimuth  and 
elevation,  g|.((l),0)  denotes  the  receive  element  antenna  power  pattern 
as  a  function  of  azimuth  and  elevation,  Nq  denotes  the  receiver 
noise  power  spectral  density,  Lg  denotes  the  system  losses,  and 

denotes  the  range  from  the  platform  to  the  clutter  ring 
illuminated  by  the  mainbeam.  Parameters  Pj,  T^,  Lg,  and  are 

radar  system  and  scenario  parameters  provided  as  inputs  to  the 
software;  the  remaining  parameters  are  calculated  as  described 
next.  Pre-specif ied  physical  constants  are  listed  in  Table  2-2. 


The  radiation  wavelength  is  determined  given  the  speed  of 
light,  V,  and  the  radiation  frequency,  fQ,  as 


For  V  in  km/sec  and  (q  in  MHz,  the  wavelength  is  calculated  in  mm. 

And  the  receiver  noise  power  spectral  density  is  defined  as 
(Skolnik,  1980) 

(2-35)  No  =  kBT^F, 

where  is  the  system  noise  figure,  a  radar  system  parameter 
provided  as  input  to  the  software.  Boltzmann's  constant,  kg,  and 
the  standard  temperature,  T^,  are  as  specified  in  Table  2-2. 
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SYMBOL 

DESCRIPTION 

VALUE 

V 

Speed  of  light 

2.998x10®  km/sec 

To 

Mean  radius  of  theEarth 

6,368.0  km 

''e 

Effective  radius  of  the  Earth  (=  (4/3)  ro) 

8,490.66667  km 

Boltzmann's  constant 

1.38x10-23  W/(Hz-K) 

To 

Standard  temperature  for  radar  system 

« 

290  K 

Yc 

Constant-gamma  clutter  reflectivity 
model  parameter  for  heavy  land  clutter 
case 

-3  dB 

Table  2-2.  Physical  constants  for  multichannel  signal  model. 


The  clutter  patch  RCS,  determined  as 

(2-36)  Cc  =  So  Ac  =  So  WcLc 

where  Sq  is  the  clutter  patch  reflectivity,  is  the  clutter  patch 
area,  is  the  clutter  patch  width,  and  is  the  clutter  patch 

length.  Clutter  patch  reflectivity  is  determined  according  to  the 
constant-gamma  model,  which  is  defined  as  (Nathanson,  1991) 

(2-37)  So  =  YcSin(¥c) 

with  the  reflectivity  model  parameter  as  specified  in  Table  2-2, 

and  the  clutter  ring  grazing  angle  \l/.  as  defined  in  Figure  2-1  and 

c 
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Table  2-1.  For  each  patch  in  the  360  deg  clutter  ring,  clutter 
patch  width  is  calculated  as 

(2-38)  Wc  =  A(l)  = 

with  the  range  to  clutter  ring  as  defined  in  Figure  2-1,  and  the 
clutter  patch  angular  extent  Acj)  as  defined  in  Table  2-1.  Also, 

is  the  number  of  clutter  patches  in  the  clutter  ring  (Table  2-1) . 
Lilcewise,  each  patch  in  the  360  deg  clutter  ring  is  of  length  Lj,, 

calculated  as 


(2-39) 


Lc  = 


V 


cos(vi/  )  2feCOs(\|/  ) 


2  tan 


rp 


1e 


sin(\j/  ) 


if  P^tan(vJ  <  rctan(|3^J 


if  P^tan(v|/J  >  r^tan(p,J 


with  p^  denoting  the  radar  range  resolution,  fg  denoting  the 
receiver  bandwidth,  p.|g  denoting  the  transmit  array  one-way,  3-dB 
elevation  beamwidth,  and  the  other  parameters  are  as  defined 
previously.  The  receiver  bandwidth  is  a  radar  system  parameter 
provided  as  input  to  the  software;  whereas  the  elevation  beamwidth 
is  approximated  as  (Skolnik,  1970) 

00  <  60° 


00  =  90° 


(2-40a) 


Ple- 


bw 


(Je-l)dcos(0o) 
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0.866  for  uniform  taper 

(2-40C)  b^=< 

1. 3  for  Dolph-Chebyshev  taper 

where  Jg  is  the  number  of  elevation  axis  elements  beamformed  into 

one  azimuth  channel,  d  is  the  linear  array  inter-element  spacing 
(assumed  to  be  the  same  for  elevation  and  azimuth)  ,  0q  is  the 

clutter  ring  elevation  angle  calculated  as  in  Table  2-2,  and  bw  is 
a  constant  factor.  Azimuth  angle  and  inter-element  spacing  are 
illustrated  in  Figure  2-1.  The  cosine  factor  in  Equation  (2-40a) 
accounts  for  scanning-induced  broadening  of  the  main  beam 
(Skolnik,  1970).  The  constant  factor  in  Equation  (2-40a)  is 
introduced  herein  to  account  for  taper-induced  beam  broadening. 
Its  minimum  value  is  0.866  for  the  uniform  taper,  as  indicated  in 
Equation  (2-40c)  .  The  value  set  in  Equation  (2-40c)  for  Dolph- 
Chebyshev  tapering  is  a  representative  value.  An  expression  for 
the  exact  value  is  unknown  since  the  beam  broadening  induced  by 
the  Dolph-Chebyshev  taper  is  a  complicated  function  of  sidelobe 
level  and  number  of  elements.  Other  elevation  pattern  tapers 
applied  to  reduce  the  antenna  sidelobes  also  result  in  a  broader 
beamwidth  (with  respect  to  the  uniform  taper) . 

The  2-D  transmit  array  has  an  antenna  power  pattern  denoted 
as  G{(<t),0).  For  simplicity,  the  elevation  dimension  is  ignored 

herein,  and  the  clutter  return  is  calculated, for  the  mainbeam 
elevation  (0q)  only.  Thus,  the  antenna  power  pattern  used  to 
generate  the  ground  clutter  model  is  Gt((|),0o).  From  standard  array 

theory,  the  antenna  pattern  is  the  product  of  the  element  pattern 
and  the  array  pattern  (Skolnik,  1980)  .  Let  denote  the 

element  power  pattern,  which  is  assumed  to  be  a  squared  cosine  in 

the  azimuth  axis, 
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cos^{(l))  for  -90°<<})<90° 

(2-41)  ge(<l).eo)  = 

G[^cos^((t))  for  90°<<t)<270° 


where  is  the  backlobe  region  attenuation.  This  element  pattern 
model  is  used  by  Ward  (1994)  .  Now  let  Ga((|),0o)  denote  the  array 

pattern  in  the  azimuth  axis,  which  is  assumed  to  be  either  uniform 
or  Dolph-Chebyshev.  Other  pattern  options  can  be  added  to  the 
model  in  the  future.  The  uniform  array  pattern  is  generated  as 
(Skolnik,  1980) 


(2-42) 


for  -90°<(|)<90° 


for  90°<  (1x270° 


where  <1)^  is  the  mainbeam  azimuth  angle  (see  Figure  2-1  and  Table  2- 
1) ;  and  the  Dolph-Chebvshev  array  pattern  (characterized  by  equi- 
level  sidelobes)  is  generated  as 

-l|F[ej(f3o)(8)u]f  for  -90°<(t)<90° 

(2-43a)  Ga((|),eo)- 

Gb-^K[§Ao)®yir  90°<())<270° 
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1 

(2-43b)  ej(f^„)=  . 

giSifsoK 

{2 -44c)  fgo  =— cos(0o)sin((t)Q) 

Aq 

here  u  is  the  J-dimensional  vector  of  Dolph-Chebyshev  weights,  F 
denotes  the  discrete  Fourier  transform  (DFT)  operator,  and  ® 
denotes  the  Schur  (element-by-element)  product.  The  Dolph- 
Chebyshev  weights  are  calculated  as  in  (Urkowitz  et  al.,  1973), 
which  requires  specification  of  the  number  of  weights  (dimension 

of  the  vector  u)  and  the  desired  sidelobe  level.  The  spatial 
frequency  vector  ej(tso)  steers  the  array  pattern  main  beam  to  the 

desired  azimuth  and  elevation,  and  the  real -valued,  scalar  U  is  a 
normalization  factor  selected  so  that  Ga((t>Q,0o)  =  1  •  Given  these 

definitions,  the  transmit  antenna  power  pattern  is  obtained  as 
(2-45)  G{{((),0q)  =  Gq  G3((t),0Q)  9e(*t*’®o^ 

where  Gq  is  the  transmit  pattern  gain.  The  gain  Gq  is  a  radar 
system  parameter  provided  as  input  to  the  software. 

As  assumed  for  the  transmit  pattern,  the  elevation  dimension 
is  ignored  herein  for  the  receive  element  antenna  power  pattern, 
gr(<t),0o),  and  the  clutter  return  is  calculated  for  the  mainbeam 
elevation  (0  =  0q)  only.  Thus,  the  receive  element  antenna  power 

pattern  is  obtained  as 
(2-46)  gr(4).eo)  =  ge(<f>-0o) 
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This  model  is  adopted  for  the  case  of  a  linear  receive  array,  as 
well  as  for  the  case  of  a  2-D  receive  array  with  beamforming  in 
the  elevation  axis  (since  the  azimuth  axis  is  unchanged) . 

The  receive  column  pattern  beamforming  gain,  G^r,  is  included 
in  the  CNR  model  to  account  for  the  cases  where  the  receive  array 
is  a  2-D  (multi-column,  multi-row)  array  of  elements  that  is 
"collapsed"  into  a  linear  (single-row)  array  by  beamforming  in  the 
elevation  (column)  axis.  Therefore,  this  gain  is  computed  as 

(2-47)  Gcr  =  JeGe 

where  Jg  is  the  number  of  elevation  axis  elements  beamformed  into 
one  azimuth  channel,  and  Gg  is  the  receive  element  gain.  Notice 
that  for  Jg  =  1  the  product  Gg^  gr(<|),6o)  =  Gg  ge(<t>t6o)  and  represents  a 
single  receive  element,  as  should  be  the  case.  Both  constants  Jg 
and  Gg  are  radar  system  parameters  provided  as  inputs  to  the 
software.  This  completes  the  CNR  model  definition. 


Consider  now  the  vector  spatio-temporal  ground  clutter 
process,  {c(n)},  as  defined  in  Equation  (2-9)  .  This  vector  process 
is  modeled  herein  as  the  sum  of  the  return  from  all  the  clutter 
patches  in  the  clutter  ring  illuminated  by  the  mainbeam.  That  is. 


Ne-1 


(2-48a)  C^+i(n)  = 

p=0 


n  =  0, 

k  =  0, 1 . J-1 


(2-48b) 


No-1 


c(n)  =  ^ 


j2Kf 


cdp 


p=0 


where  denotes  the  complex-valued  clutter  return  from  the  single 
scatterer  representing  the  pth  patch,  and  f^dp  and  denote  the 
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normalized  Doppler  frequency  and  spatial  frequency,  respectively, 
of  the  pth  clutter  patch.  These  frequencies  are  computed  as 

(2-49)  fcdp  =  0.5  cos(eQ)  sin((t)p  -  y) 

(2-50)  fcsp  =  (ci/X(,)  005(0^)  sin((l)p)  =  0.5  cos(0o)  sin((l)p) 

(2-51)  (t)p  =  -  71  +  p  A(t)  p  =  0, 1, . . . ,  Nc- 1 

with  V  the  clutter  ridge  slope  as  defined  in  Equation  (2-32); 
also,  Acj)  is  the  clutter  patch  angular  extent,  and  N^.  is  the  number 

of  clutter  patches  in  a  clutter  ring  (see  Figure  2-1  and  Table  2- 
1) .  Nc  is  a  scenario  parameter,  and  is  provided  as  input  to  the 

software.  Notice  that  the  clutter  patch  spatial  frequencies  are 
bounded  by  ±0.5;  in  contrast,  the  clutter  patch  Doppler  frequencies 
are  unbounded.  As  discussed  earlier,  Doppler  frequencies  greater 
than  0.5  or  smaller  than  -0.5  are  ambiguous,  and  these  frequencies 
fold  over  onto  the  observable  frequency  range,  -0.5<fd<0.5. 

The  clutter  return  from  the  single  scatterer  representing  the 
pth  patch  is  modeled  herein  as  a  complex-valued,  scalar  random 
variable  with  mean  zero  and  Gaussian— distributed  real  and 
imaginary  parts  (Rayleigh-distributed  magnitude) .  Also,  all 
clutter  patches  are  assumed  to  be  pairwise  statistically 
independent.  It  follows  from  the  definition  of  CNR  that  the  power 
(variance)  of  the  pth  patch  clutter  return,  Ocp '  obtained  as 


(2-52) 


ElVqj 


4Vq='’n  y+p.eolV., 


-71  <  <t)p  <  1 .571 


where  6p.q  is  the  Kronecker  delta,  which  is  defined  as 
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1  if  p  =  q 

(2-53)  Vq=‘ 

0  otherwise 


The  receiver  noise  variance,  ,  is  a  radar  system  parameter,  and 

is  provided  as  input  to  the  software.  With  these  definitions,  the 
clutter  true  auto-covariance  matrix  sequence  {R  CC  (m)}  is  obtained  as 


(2-54)  Rcc('T^)  =  E[c(n)cH(n-m)]  = 


^cp  ^ 


-j(^csp  )  §J  (^csp ) 


m>0 


with  negative  lags  obtained  as 

(2-55)  Rcc(-m)  =  E[c(n)c^^(n+m)]  =  Rcc(m)  m>0 

Each  matrix  in  the  sequence  is  JxJ.  Matrix  Rcc(fTi)  is  Hermitian  only 
for  m  =  0;  however,  matrix  Rcc(m)  has  Toeplitz  structure  for  all  m 
due  to  the  spatially-induced  relationship  among  the  elements  of 
c(n),  as  expressed  in  Equation  (2-9) .  Equation  (2-54)  is  analogous 
to  Equation  (2-24) ,  and  the  procedure  used  to  generate  the  signal 
covariance  matrix  sequence  can  be  used  to  generate  the  clutter 
covariance  matrix  sequence. 


2 . 6  riiannel  Output  True  Covariance  Model 


The  channel  output  vector  sequence  {x(n)}  is  modeled  as  the  Siam 
of  the  signal  (moving  target),  broadband  interference  (jamming), 
ground  clutter,  and  receiver  array  noise  (Equation  (2-3)).  These 
four  components  are  assumed  to  be  pairwise  statistically 
independent.  Thus,  the  channel  output  true  covariance  matrix 
sequence  is  obtained  as 

(2-56)  Rxx(fn)  =  E[x(n)  x^(n-m)  ]  =  Rss(m)  +  Rii(m)  +  Rcc(m)  +  Rww(fn)  ^  0 
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with  negative  lags  obtained  as 


(2-57)  Rxx(-m)  =  E[x(n)x^(n+m)]  =  R^^(m)  m>0 

Matrix  Rxx(*^)  Toeplitz  for  all  values  of  m  because  each 

individual  matrix  on  the  right-hand-side  of  Equation  (2-56)  is 
Toeplitz  or  zero  for  all  ID  >  0.  Likewise,  matrix  Rxx(O)  is 

Hermitian  because  each  individual  matrix  on  the  right -hand- side  of 
Equation  (2-56)  is  Hermitian  for  171  =  0. 

2 . 7  Signal  (Movincr  Target)  Data  Generation 

Moving  targets  are  modeled  as  defined  in  Equations  (2-15) -(2- 
21)  of  Section  2.1.2.  That  is,  as  the  output  vector  sequence  of  a 
linear,  stationary,  discrete-time,  complex-valued  system  in  state- 
space  form  driven  by  a  complex-valued,  zero-mean,  Gaussian- 
distributed  white  noise  vector  sequence,  {u.(n)},  with  diagonal 

covariance  matrix  Q.  The  initial  condition  for  Equation  (2-15a) , 
y(0),  is  a  complex-valued,  zero-mean,  Gaussian-distributed  Gg- 

dimensional  vector.  This  model  is  a  set  of  Gg  scalar  systems  in 

parallel,  and  each  scalar  system  represents  a  moving  target.  Data 
generation  is  accomplished  by  driving  the  system  (2-15a)  with  a 
pseudo-random  initial  condition,  and  a  pseudo-random  input 
sequence  of  N-point  duration.  The  initial  condition  is  a  random 
sample  from  a  Gaussian  distribution  with  mean  vector  zero  and 
covariance  matrix  P.  Matrix  P  is  the  steady-state  covariance 
matrix  of  the  state  vector,  and  it  satisfies  a  Lyapunov  equation 
of  the  form 

(2-58)  P  =  FPF^  +  Q 

For  an  asymptotically-stable  system,  Equation  (2-58)  has  a 
positive-definite  solution.  Furthermore,  since  matrices  F  and  Q 
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are  both  diagonal,  Equation  (2-58)  reduces  to  the  following  set  of 
scalar  equations, 

(2-59)  P(k,k)  =  IFP  P  +  Q(k,k)  k=1,2 . Gg 

which  is  an  equivalent  expression  for  Equation  (2-19b) .  Selection 
of  P  as  the  covariance  matrix  of  the  state  vector  initial 
condition  places  the  system  model  in  stochastic  steady-state  at 
time  n  =  0,  and  it  remains  in  stochastic  steady-state.  The  kth 
moving  target  contributes  power  to  each  channel  output,  and  the 

total  signal  power  in  each  channel  is 

(2-60)  af  =  diag(P)  = 

k=1 


with  P(k,k)=a||^,  as  defined  in  Equation  (2-21). 

Swerling  Cases  1  and  2  can  be  generated  by  appropriate  choice 
of  the  model  parameters  in  model  (2-16) - (2-18) .  For  notational 
simplicity,  consider  the  case  of  a  single  target,  k  =  1.  For 
Swerling  Case  1,  let  =  1  and  ^d1  “  which  forces  the  driving 

noise  covariance  to  be  equal  to  zero,  and  the  state-space  model 
becomes 

Swerling  Case  1: 


(2-61a) 

y(n+1)  =y(n) 

0<n<N- 1 

(2-61b) 

s(n)  =  ej(f,3i)  y(n) 

0  <  n  <  N  - 1 

(2-61c) 

y(0)  ~  mo, 0.5  4 12) 
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where  0.5  Ogjg)  denotes  the  2-D  circular  Gaussian  distribution 

for  a  complex-valued  scalar  with  mean  zero  and  diagonal  covariance 
(the  scalar,  complex- valued  initial  condition  has  variance  . 

For  Swerling  Case  2,  let  =  0,  which  forces  the  driving  noise 

covariance  to  be  equal  to  unity,  and  the  state-space  model  becomes 


Swerlina 

Case  2 : 

(2-62a) 

y(n+1)  =  u(n) 

0  <  n  <  N  - 1 

(2-62b) 

s(n)  =  ej(ft3i)y(n) 

0  <  n  <  N  - 1 

(2-62c) 

u(n)  ~  fA^(0,0.5l2) 

0  <  n  <  N  - 1 

(2-62d) 

y(0)  ~  mms  I2) 

Swerling  Cases  3  and  4  can  be  modeled  with  the  equations  for  Cases 
1  and  2,  respectively,  if  the  initial  condition  vector  and  the 
input  noise  vector  are  distributed  as  chi-square  with  four 
degrees-of- freedom  instead  of  Gaussian. 

A  note  of  caution  is  appropriate.  Swerling  Case  1  (also  Case 
3)  conditions  correspond  to  a  system  model  which  is  conditionally 
stable,  as  opposed  to  asymptotically  stable.  This  violates  one  of 
the  conditions  for  applicability  of  the  identification  algorithms 
used  in  this  program  (as  well  as  many  other  identification 

algorithms),  and  can  lead  to  difficulties.  Specifically,  it  has 
been  observed  in  various  simulations  that  setting  a|^  =  1  leads  to 

canonical  correlations  greater  than  unity,  which  is  an  invalid 
condition.  The  exact  theoretical  answer  for  such  cases  is  unity 
canonical  correlations,  but  numerical  issues  often  force  the 
invalid  result. 
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2 . 8  Broadband  Interference  (Jamming)  Data  Generation 


Broadband  interference  (jamming)  is  modeled  as  defined  in 
Section  2.1.3.  That  is,  the  broadband  interference  sequence,  {i(n)}, 
is  a  complex-valued,  Gaussian-distributed,  zero-mean,  white 
sequence,  which  represents  the  received  interference  from  all 
direct,  barrage,  point-souce  jammers.  In  this  model  each  channel 
receives  the  same  amount  of  power  from  each  jammer.  It  follows 
that  the  jamming  sequence  model  is 

Gi 

(2-63a)  i(n)  =  2^  ik(n)  ej(fj3k)  0<n<N-1 

k=1 

(2-63b)  ik(n)  ~  (0,0.5 afjg)  1  <  k  <  Gji  0  <  n  <  N  - 1 


where  ik(n)  is  the  kth  jamming  source  white  sequence  with 
corresponding  spatial  frequency  fjgk,  and  is  the  variance  (power) 

of  the  kth  jamming  source  at  each  channel  output.  Data  generation 
is  accomplished  by  generating  G|  complex-valued,  pseudo -random, 
scalar,  white  sequences  {ik(n)}  according  to  Equation  (2-63b)  ,  and 

combining  them  according  to  Equation  (2-63a) .  The  total  jammer 
power  in  each  channel,  denoted  as  Oj  ,  is  given  as 


(2-64) 


for  Gj  independent  interference  sources . 

2 . 9  Receiver  Array  Noise  Data  Generation 

Receiver  array  noise  is  modeled  as  described  in  Section 
2.1.4.  That  is,  the  receiver  array  noise  sequence,  {w(n)},  is  a 
complex-valued,  Gaussian-distributed,  zero-mean,  white  sequence 
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uncorrelated  from  channel  to  channel  (with  diagonal  covariance  at 
each  time  instant) , 


(2-65a) 

E[w(n)w>-m)]  = 

(2-65b) 

w(n)  ~  fA^(0.0.5a2y 

0  <  n  <  N  - 1 

where 

is  the  noise  power  in  each  channel. 

Data  generation  is 

accomplished  by  generating  a  J-dimensional,  complex- valued,  pseudo¬ 
random,  white  vector  sequence  according  to  Equation  (2-65) . 

2.10  Ground  Clutter  Data  Generation 

Ground  clutter  is  modeled  as  described  in  Section  2.1.5. 
That  is,  the  ground  clutter  sequence,  {c(n)},  is  the  sum  of  the 
return  from  ground  patches  in  the  clutter  ring  illuminated  by 

the  array  main  beam.  Each  patch  is  represented  by  a  single 
complex-valued,  Gaussian-distributed,  zero-mean  random  scatterer, 
as  indicated  in  Equation  (2 -48b) ,  and  the  scatterers  are  assumed 
to  be  pairwise  statistically-independent .  Equation  (2-48b)  is 
repeated  next  as  Equation  (2-66) ,  for  convenience, 

Nc-1  . 

(2-66)  SW  =  "  =  . N-1 

p=0 

(2-67)  rip  ~  5\/:(0,0.5a2p  Ig)  p  =  0,1 . 

As  before,  rip  denotes  the  clutter  return  from  the  single  scatterer 
representing  the  pth  patch,  and  and  denote  the  normalized 
Doppler  frequency  and  spatial  frequency,  respectively,  of  the  pth 
clutter  patch.  These  frequencies  are  computed  as  in  Equations  (2- 
49) -(2-50).  The  clutter  process  {c(n)}  has  true  covariance  matrix 
sequence  {Rcc(m)}  defined  as  in  Equations  (2-54)  and  (2-55) .  Data 
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generation  is  accomplished  by  generating  complex-valued, 

pseudo -random,  scalar,  random  variables  {T|p}  according  to  Equation 

(2-67),  and  combining  them  according  to  Equation  (2-66).  This 
procedure  is  computationally  intensive  when  the  number  of  ground 
patches  (N(.)  is  large  and  the  number  of  receive  array  elements  (J) 
is  large.  The  pth  scatterer  contributes  power  to  each  channel 

output,  and  the  total  ground  clutter  power  in  each  channel  is 
Nc 

(2-68)  a2=^afp 

p=i 


Equation  (2-68)  follows  from  Equations  (2-48a)  and  (2-52)  . 

2.11  Channel  Output  Data  Generation 

The  channel  output  vector  sequence  {x(n)}  is  modeled  as  the  sum 
of  the  signal  (moving  target) ,  broadband  interference  (jamming) , 
ground  clutter,  and  receiver  array  noise.  These  four  components 
are  assumed  to  be  pairwise  statistically  independent.  Thus,  the 
channel  output  vector  sequence  {x(n)}  is  generated  as  the  sum  of  the 
four  independent  components, 

(2-69)  x(n)  =  s(n)  +  i(n)  +  c(n)  +  w(n)  0<n<N-1 

and  the  total  power  in  each  channel  output  is 

O  p  2  p  p 

(2-70)  =  CTg  +  CTj  +  (Jq  +  Gp 

Single-channel  power  ratios  (such  as  signal-to-noise  ratio)  can  be 
calculated  given  the  channel  output  power  level  for  each 
individual  component  (signal;  interference;  clutter) . 
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2.12  Covariance.  Spectrum,  and  Graphical  Issues 

Diagnostics  constitute  a  key  part  of  the  software  model 
developed  in  this  program.  Generation  of  the  2-D  power  spectrum 
of  the  surveillance  radar  space/ time  process  and  plots  of  the 
spectrum  and  associated  auto-covariance  sequence  provide  insight 
into  the  structure  and  characteristics  of  the  "true"  process  and 
the  simulated  data.  Power  spectrum  and  auto-covariance  matrix 
estimates  are  important  also.  Those  issues  are  discussed  in  this 
section. 

The  channel  output  vector  sequence  {x(n)}  can  be  viewed  as  a 
scalar  2-D  sequence  due  to  the  relationship  that  exists  between 
elements  of  the  vector  processes  s(n),  c(n),  and  i(n),  and  because  the 
spatially-uncorrelated  white  noise  vector  w(n)  can  be  viewed  also 
as  a  scalar  2-D  white  sequence.  Specifically, 

(2-71)  x,^^^(n)  =  x^(n+kx)  =  x.|(n,k)  0  <  n  <  N-1;  0  <  k  <  J  - 1 

The  scalar  2-D  process  {x^(n,k)}  has  a  scalar  2-D  power  spectrum 
denoted  as  {Sxx(fd>y}'  ^  scalar  2-D  auto-covariance  sequence 
(ACS)  denoted  as  {rxx(m,^)},  where  (m,^)  denotes  the  lag  pair 

corresponding  to  the  frequency  pair  (fjj.fg).  The  2-D  power  spectrum 
and  the  ACS  form  a  Fourier  transform  pair  (see,  for  example. 
Dudgeon  and  Mersereau,  1984), 


oo  oo 

(2-72)  Sxx(fd’U  =  ^['■xx(m.^)]  =  exp[-j27t(fdm  +  fg^)] 

m=-«>  f=-oo 


-  0.5  <  <  0.5  and  -  0.5  <  fg  <  0.5 
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(2-73) 


hxM)  =  ^'^[Sxx(fd>y]  = 


// 


Sxx(fd-^s)  exp[j27c(fdm  +  fg^)]  df^  dfg 


-  <  m  <  03  and  -  oD  <  /  <  oo 


In  this  transform  pair  the  ACS  is  infinite  and  the  power  spectrum 
is  continuous.  Neither  of  these  two  conditions  is  true  in 
practical  applications  such  as  surveillance  radar.  In  fact,  the 
number  of  lags  available  is  determined  by  the  number  of  channels 
and  the  number  of  pulses  in  one  CPI,  and  the  spectrum  is  sampled 
at  sampling  intervals  determined  also  by  the  number  of  channels 
and  the  number  of  pulses  in  one  CPI.  For  these  conditions  the 
relation  that  applies  is  the  DFT  and  its  inverse;  namely. 


N-1  J-1 


(2-74a) 

Sxx(fd-fs)  =  F[rxx(m.-^)]  = 

X  X  ®xp[- j27t(fdm  -f  fg  £)] 

m=-N+1  ^=-J+1 

(2-74b) 

fd  =  kd^^d 

-  N  +  1  <  k J  <  N  - 1 

(2-74C) 

fs  =  ksAfs 

-  J  1  <  kg  <  J  - 1 

(2-74d) 

“  2N  - 1 

{2-74e) 

"  2J-1 

{2-75a) 

r,x(nn,fl  =  F-'(S,,(f,,g] 

-  N -1- 1  <  m  <  N  - 1; 

0.5  0.5 

(2-75b) 

rxx(m.0=  X  X  SxxMs)®>^P027t(fdfTi+y)]  AfdAfg 

f  =-0.5  f,=-0.5 
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0.5  0.5 

(2-750  =  — -L—  X  S  S,,(V3)exp[j27c(fdm  +  f3^)] 

'  '  f^=-0.5  f^=-0.5 

Parameters  Af^  and  Afg  determine  the  spectrum's  frequency  resolution 

along  the  Doppler  and  spatial  axes,  respectively.  Periodic 
sampling  of  the  continuous  spectrum  causes  the  ACS  to  be  repeated 
periodically  in  the  lag  domain.  Likewise,  the  sampled  spectrum  is 
repeated  periodically  in  the  frequency  domain  because  the  ACS  is 
available  only  at  discrete  times . 

It  is  possible  to  express  the  2-D  ACS  in  terms  of  the  2-D 
sequence  {x^(n,k)}  as  follows: 

(2-76)  rxx(m,^)  =  E[x/n,k)Xi(n-m,k-^)]  m>0:  ^>0 

with  negative  lag  values  determined  according  to 
(2-77)  =  ^>0;  ^>0 

(2-78)  =  m>0:  ^>0 

(2-79)  rxx(-m,-^)  =  r*x(m,^)  m>0:  ^>0 

The  channel  output  true  covariance  matrix  sequence,  {Rxx(m)}»  can  be 
expressed  also  in  terms  of  the  channel  output  2-D  ACS  using 
Equations  (2-76) - (2-79)  .  The  JxJ  covariance  matrix  sequence  for 
the  process  {x(n)}  is  defined  in  Equations  (2-54)  and  (2-55)  .  For 
each  m,  matrix  Rxx(f^)  can  be  expressed  in  terms  of  the  elements  of 

the  2-D  ACS  as 
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(2-80) 


Rxx(m)  = 


rxxlm.O)  rxx(m,-1) 


rxx(m.1) 


rxx(m,0) 


rxx(m,J-2)  rxx(m,J-3) 
rxx(m,J-1)  rxx(m.J-2) 


rxx(m,-J  +  2)  rxx(m,-J  +  1) 
rxx(m,-J+3)  rxx(m,-J+2) 


rxx(m.O)  rxx(m,-1) 

TxxCm.l)  rxx(m,0) 


Notice  that  Rxx(m)  is  defined  completely  by  the  first  row  and  the 

first  column,  as  befits  its  Toeplitz  structure.  Based  on 
Equations  (2-77 )- (2-79 )  ,  Rxx(hl)  is  Hermitian  only  for  m  =  0,  as 

stated  previously. 


Equation  (2 -74a)  can  be  expressed  as  a  matrix  equation  of  the 

form 

(2-81)  Sxx  =  TlRxxTr 

where  S^^x  i®  (2J-1)x(2N-1)  power  spectrum  matrix,  Tl  is  the  (2J- 
1)x(2J-1)  DFT  matrix,  Tr  is  the  (2N-1)x(2N-1)  DFT  matrix,  and  Rxx  is 
the  (2J-1)x(2N-1)  ACS  matrix  defined  as 


(2-82) 


-rxx(-N  +  1 -J  +  1)  - 

*^xx(  J  +  1) 

^xx(®>  J  +  1) 

«'xx(N-1,-J  +  1) 

fxx(“N  +  1,-1) 

*xx(®>”l) 

rxx(N-1,-1) 

r„(-N  +  1,0)  ■■■ 

rxx(-1,0)  ^ 

^xx(^>®) 

r„(N-1,0) 

fxx(“N  +  1,  J-1) 

rxx(0,J-1)  - 

r„(N-1,J-1) 

DFT  matrices  T|_  and  Tr  operate  on  the  columns  and  rows, 
respectively,  of  the  ACS  matrix  Rxx  to  implement  the  2-D  DFT.  It 

is  assumed  herein  that  the  DFT  matrices  also  include  any  row 
and/or  column  shift  operations  necessary  to  generate  a  spectrum 
matrix  with  origin  at  the  (J,N)th  element.  Additionally,  it  may  be 
necessary  to  take  the  absolute  value  of  the  right-hand-side 
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Equation  (2-81)  in  order  to  obtain  the  true  real-valued  spectrum; 

this  is  a  function  of  the  way  a  specific  DFT  is  implemented.  The 
partitions  on  in  Equation  (2-82)  represent  the  four  quadrants, 

Q1  through  Q4,  of  the  2-D  plane,  with  the  lower-right  partition 
corresponding  to  Ql;  specifically. 


(2-83) 


'XX 


Q2  !  Q4 
Qs]  Ql 


Notice  that  the  origin  is  assigned  to  Ql .  This  agrees  with  the 
convention  adopted  in  most  three-dimensional  (3-D)  software 
graphics  paclcages,  including  MATLAB.  The  quadrants  are  shown  in 
Figure  2-3,  where  the  axes  are  drawn  in  agreement  with  MATLAB 's  3- 
D  graphics  default  axes  orientation  and  viewing  angles. 


Figure  2-3.  MATLAB  two-dimensional  grid  definition  for  default 
view  angles  (azimuth  and  elevation) . 
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In  the  approach  formulated  herein  for  the  channel  output  true 
covariance  model,  the  covariance  matrix  sequence  is  generated 
directly,  without  having  to  generate  the  power  spectrum.  This  is 
appropriate  for  generating  the  inputs  required  by  identification 
and  detection  routines  when  running  "known  covariance"  cases. 
However,  it  is  difficult  to  appreciate  the  features  of  the  channel 
output  by  examination  of  either  the  true  ACS  or  the  true 
covariance  matrix  sequence.  In  contrast,  visual  examination  of  a 
3-D  plot  of  the  power  spectrum  is  very  informative. 

Given  a  true  ACS  or  covariance  matrix  sequence,  the  (2J-1)x(2N- 
1)  power  spectrum  matrix  is  obtained  via  Equation  (2-81)  ,  and  a 

3-D  plot  can  be  generated  at  frequency  resolutions  determined  by 
Equations  (2-74d)  and  (2-74e) .  The  resulting  plot  will  exhibit 
the  principal  features  present  in  the  true  covariance  model,  but 
will  have  three  drawbacks.  First,  the  frequency  resolution 
specified  by  Equations  (2-74d)  and  (2-74e)  is  the  coarsest 
resolution  possible  (determined  by  the  fixed  values  J  and  N)  . 
Second,  the  frequency  leakage  level  resulting  from  the  Fourier 
transformation  in  Equation  (2-81)  is  the  highest  possible  because 
the  true  ACS  lags  beyond  those  specified  in  Equation  (2-75)  are 
unavailable,  which  is  equivalent  to  multiplying  the  infinite- 
duration  true  ACS  by  a  rectangular  lag  window  defined  with  its 
corners  at  the  lag  pairs:  (N-1,J-1),  (-N+1,J-1),  (-N+1,-J+1),  and  (N-1,- 
J+1).  Third,  for  large  values  of  N  and/or  J,  the  numerical 
evaluation  of  Equation  (2-81)  is  inefficient  because  the  length  of 
the  required  transforms  are  odd  numbers  (2N-1  and  2J-1  )  ,  and 
efficient  DFT  algorithms  are  unavailable,  in  general,  for  large 
and  odd  length  values . 

Each  of  the  drawbacks  mentioned  above  can  be  addressed,  as 
shown  next.  The  first  drawback  is  overcome  by  "padding"  the  true 
ACS  with  zeros .  The  second  drawback  is  mitigated  by  applying  a 
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low-sidelobe  lag  window  to  the  finite-duration  true  ACS.  Use  of  a 
lag  window  improves  frequency  leakage  but  introduces  a  loss  in 
frequency  resolution  (with  the  result  that  features  extend  over  a 
larger  frequency  band) .  Finally,  the  third  drawback  is  eliminated 
by  forcing  the  length  of  the  zero-padded  ACS  along  each  axis  to  be 
an  even  integer  that  is  a  power  of  two,  which  allows  use  of  Fast 

Fourier  Transform  (FFT)  algorithms.  To  express  these  steps 
mathematically,  first  define  an  NjxN,  zero-padded  and  windowed  ACS 

matrix  as 


(2-84) 


W0R^^  [0] 

[0]  [0] 


where  W  is  a  (2J-1)x(2N-1)  lag  window  matrix  and  0  denotes  the  Schur 
product,  as  before.  The  null  sub-matrices  in  Equation  (2-84), 
denoted  as  [0]  ,  represent  appropriately-dimensioned  matrices  of 
zeros.  The  integer  Nj  is  selected  to  be  a  power  of  two  such  that 
N,  >  max(2J-1,2N-1).  Then  define  an  N,xNf  high-resolution  power 

spectrum  matrix  as 


(2-85) 


^xx  - 


\*^xx'^R 


DFT  matrices  T|_  and  Tp  in  Equation  (2-72)  are  both  NjxNj,  and  as 

before,  include  any  row  and/or  column  shift  operations  necessary 
to  generate  a  spectrum  matrix  with  origin  at  the  (N/2,N/2)th 

element.  The  absolute  value  operation  in  Equation  (2-85)  insures 
that  the  spectrum  is  real-valued,  and  is  necessary  in  this  case 
because  the  enhanced  ACS  matrix,  R^,  lacks  symmetry  with  respect 

to  its  middle  element,  which  generates  complex-valued  terms  in  the 
product  T|_Rj(^Tp  (this  effect  is  common  to  most  DFT  implementations, 

and  taking  the  absolute  value  generates  the  desired  result) .  The 
frequency  resolution  associated  with  spectrum  Sj^  is 
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(2-86) 


AL  =  Af  =  Af  =  — 

®  N, 

Finally,  the  lag  window  matrix  W  is  selected  as  the  (2J-1)x(2N-1) 

matrix  formed  from  a  real-valued,  finite-duration,  separable,  2-D 
lag  window  sequence  {a)(nn,^)};  that  is, 

(2-87)  (o(m,^)  =  a)i(m)  coi(^)  -  N  +  1  <  m  <  N  - 1;  -  J  +  1  <^ <  J  - 1 


^(-■N  + 1,  —  J  + 1) 

•  ®(-1,-J  +  1) 

(o(0,-J  +  1)  ••• 

co(N-1,-J  +  1)' 

co(-N  +  1-1)  • 

•  co(-1-1) 

co(0-1)  - 

a)(N-1,-1) 

o)(-N  +  1,0) 

•  a(-x~0)  " 

co(0,0) 

co(N-1,0) 

co(-N  +  1,J-1)  •• 

■  (0(-1,J-1) 

co(0,J-1)  ••• 

co(N-1,J-1) 

where  {co^(*)}  is  a  real-valued,  one-dimensional  (1-D)  lag  window 

sequence.  A  1-D  lag  window  normalized  to  have  unity  value  at  the 
origin  (m  = 0)  is  preferred  because  it  preserves  total  power. 
Since  the  2-D  window  sequence  is  separable,  the  spectrum  of  the 
lag  window  matrix  (2-88)  is  circularly  symmetric. 


All  the  relations  introduced  heretofore  are  based  on  true 
quantities,  as  opposed  to  estimated  quantities,  because  only  the 
true  model  has  been  considered  up  to  this  point.  However,  several 
relations  involving  estimates  of  covariance  and  spectra  are 
relevant  for  data  generation  and  diagnostics,  and  are  introduced 

below.  Prior  to  that  it  is  convenient  to  introduce  the  concept  of 
an  extended  2-D  sequence  {x^(n,k)}  obtained  by  padding  the  sequence 
{x.|(n,k)  I  n  =  0,  1 ,  .  .  .  ,  N-1 ;  k  =  0,  1 . J-1}  with  zeros  as 


(2-89) 


Xi(n,k) 


0<n<N-1:  0<k<J-1 
N<n<2N-2;  J<k<2J-2 
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This  extended  sequence  is  important  in  ACS  and  spectrum  estimator 
definitions  established  next. 

Given  the  2-D  sequence  {x.|(n,k)},  the  non-circular,  biased, 
time-average  estimate  of  the  ACS  for  {x^(n,k)}  is  defined  as 


0  <  m  <  N  - 1 

where  the  caret  ('^)  over  a  variable  denotes  an  estimator,  or 
estimate,  of  the  variable.  A  similar  formula  (with  different 
limits  for  the  summations)  is  available  to  estimate  the  ACS  for 
negative  lags,  but  it  is  simpler  to  obtain  the  ACS  estimate  for 
negative  lags  according  to  the  equivalence  relations 


(2-91) 

0  <  m  <  N  - 1 : 

0<^<  J-  1 

(2-92) 

E 

*  X 

II 

1 

E 

X 

0  <  m  <  N  - 1; 

0<^< J - 1 

(2-93) 

r^^(-m,-0  =  rx*x(m,^) 

0  <  m  <  N  - 1; 

0<.^<J-  1 

Notice 

that  Equations  (2-91) -  (2-93 ) 

are  analogous  to 

the  negative 

lag  equations  for  the  true  ACS,  Equations  (2-77) '- (2-79)  . 

Alternative  ACS  estimators  are  obtained  via  modifications  to 
Equation  (2-90)  .  An  unbiased,  non-circular  estimator  is  obtained 
by  replacing  the  factor  JN  in  the  denominator  of  Equation  (2-90) 
by  the  factor  (J-^)(N-m).  In  the  limit  as  both  J  — >  0°  and  N— ><»,  the 

estimation  error  approaches  zero  in  an  unbiased  estimate,  but 
approaches  a  fixed,  non-zero  value  in  a  biased  estimate.  A  non¬ 
circular,  unbiased,  time-average  ACS  estimate  is  transformed  into 


N-1  J-1 


(2-90)  =  ^  X  Xi(n.k)  x*(n-m,k-.^) 


n=m  k= 
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a  biased  ACS  estimate  by  appropriately  scaling  the  unbiased 
estimate  as 

(2-94)  =  -  N  +  1  < m  <  N  - 1;  -  J  +  1  <^ <  J  - 1 

where  the  U  and  b  subscripts  denote  unbiased  and  biased  estimates, 
respectively,  since  both  types  of  estimates  appear  in  the  same 
equation.  In  subsequent  usage,  however,  the  estimator  type 
subscript  is  dropped  when  the  type  of  estimate  (or  estimator)  used 

is  clear  from  the  context  and  confusion  is  unlikely  to  occur.  The 
2-D  sequence  {(Oj(m,^)}  is  a  real-valued,  separable,  lag  window 

sequence  computed  as  the  product  of  two  normalized  1-D  triangular 
sequences, 

(2-95)  C0t(m,£)  =  coy(m)  C0f5(^)  -  N  +  1  <m  <  N  - 1;  -  J  +  1  <  J  - 1 

(2-96)  w,j(m)  =^{1,2,  ....  N-1,  N,  N-1,  ...  ,2,  1} 

(2-97)  O)t3(.0=j{1.2 . . 2,1} 

Notice  that  the  central  element  (corresponding  to  m  =  ^  =  0)  of  all 

three  sequences  is  unity.  As  a  consequence,  the  variance  (total 
power)  is  unchanged  after  the  scaling  by  {(Oj(m,^)}.  A  (2J-1)x(2N-1) 

lag  window  matrix  Wj  is  associated  with  the  2-p  sequence  {(jL)j(m,^)}, 
as  in  Equations  (2-87)  and  (2-88) . 

An  unbiased,  circular  estimator  is  obtained  by  repeating  the 
finite  sequence  over  the  time  and  space  domain,  which  results  in 
that  for  an  (N,  J)-point  sequence  all  lags  are  computed  as  the  sum 
of  JN  non-zero  products.  Of  more  relevance  to  the  approach 
pursued  herein,  Estimator  (2-90)  is  equivalent  to  a  circular , 
biased,  time-averaae  estimate  of  the  ACS  for  {x,|(n,k)}  computed  as 
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(2-98) 


0  <  m  <  N  - 1 


2N-2  2J-2 

rxx(m,^)=jL-  X  X  ><i(n>k))^(n-m,k-^) 


n=m  k=^ 


with  negative  lags  computed  as  in  Equations  (2-91) - (2-93) .  Even 
though  it  is  calculated  as  a  circular  estimator,  estimator  (2-98) 
is  equivalent  to  estimator  (2-90)  because  the  extended  sequence 
can  be  repeated  over  the  time  and  space  domain  without  altering 
the  sum  of  products,  in  relation  to  Equation  (2-90),  due  to  the 
zeros  present  in  the  extended  sequence. 

Other  ACS  estimators  can  be  defined,  and  each  one  has 
distinct  properties.  For  model  identification  estimator  (2—90)  is 
preferred  because  it  leads  to  improved  model  parameter  estimates. 
In  summary,  for  a  2-D  sequence  of  duration  (N, J),  an  ACS  estimate 
of  maximum  duration  (2N-1,2J-1)  is  obtained. 

The  2-D  neriodoaram.  denoted  as  (fj , fg ) ,  is  an  estimator  of 
the  power  spectrum  of  {x.|(n,k)},  that  is  obtained  as  the  2-D  DFT  of 
the  non-circular,  unbiased,  time-average  ACS  estimate  (see,  for 
example,  Oppenheim  and  Schafer,  1975), 


(2-99a) 

N-1 

^xxC^d’^s )  “  ~  X, 

m=-N+1# 

J-1 

X  rxx(m.^)exp[-j27t(fdm  +  fs^)] 

=-j+i 

(2-99b) 

-  N  +  1  <  kd  <  N  - 1 

(2-99c) 

fs  =  ksAfs 

-  J  +  1  <  kg  <  J  - 1 

(2-99d) 

Afd=  ^ 

^  2N-1 
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2-99e)  Af  = - 

®  2J-1 


It  follows  that  the  inverse  transform  relation  is 


2-lOOa)  rxx(m,^)  =  F'’[Sxx(fd4) 


-N  +  1<m<N-1:  -J  +  1<^<J-1 


0.5  0.5 


(2 -10 Ob)  rxx(m,.^)  =  21  Sxx(fd.fs)®xp[j27t(fdfn+fs^)]  AfjjAfj 


f^=-0.5  f^=-0.5 


0.5  0.5 


(2-lOOc)  rxx(m,-g)  =  ^2j.i)|2N-1)  ^  ^  Sxx(fd.fs)  ®xp[j27t(fjjm  +  fg^)] 


f,=-0.5  f,=-0.5 


Notice  that  transform  equation  pair  (2-99)  and  (2-100)  for  the 
estimators  is  analogous  to  transform  equation  pair  (2-74)  and  (2- 
75),  which  is  valid  for  the  true  functions. 

In  most  texts  the  periodogram  is  defined  first  as  a  function 
of  the  data  sequence,  and  then  the  equivalence  with  Equation  (2- 
99)  is  stated.  Specifically,  the  periodogram  is  computed  directly 
using  from  the  extended  2-D  sequence  as  (see,  for  example, 
Oppenheim  and  Schafer,  1975) 


(2  101a)  ^2J-1)(2N-1) 


2N-2  2J-2 


(2-lolb)  Sxx(fd.fs)  =  ^^j7^fy^^j^  S  Xxi(n.k)exp[-j27t(fdn  +  f3k)] 


(2-lOlc)  fd  =  k^Af^ 


N  +  1  <  kjj  <  N  - 1 


2-lOld)  f3  =  k3Af3 


-J  +  1  <kg<J-1 
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(2-lOle) 


1 

2N-1 


(2-lOlf) 


1 

2J-1 


The  frequency-domain  pair  and  resolution  parameters  are  repeated 
in  this  equation  in  order  to  emphasize  the  equivalence  with 
Equation  (2-99)  .  In  the  frequency  domain,  total  power  (variance) 
is  estimated  according  to  the  relation 


(2-102) 


0.5  0.5 


0.5  0.5 


^xxCO-O)  -  “  /2J.1U2N-1)  ^  ^  Sxx(Vs) 

f^=-0.5  f,=-0.5  ^  ''  '  f-=-0.5  t=-0.5 


s  'd 


Equation  (2-102)  corresponds  to  Equation  (2-100)  with  m  =  /  =  0,  and 
can  be  used  to  provide  one  indication  of  validity  of  a  spectriim 
estimate . 


The  relations  in  Equations  (2-89) - (2-101)  can  be  summarized 
in  a  diamond-shaped  diagram,  as  in  Figure  2-4.  This  diagram  is 
useful  for  visualization  of  the  analytical  issues  that  relate  the 
2-D  sequences.  The  corners  of  the  diamond,  labeled  1  through  4, 
represent  the  four  2-D  sequences,  and  the  sides  of  the  diamond 
represent  the  analytical  expressions  that  relate  the  sequences  to 
each  other.  Notice  that  only  two  operations  are  reversible,  1  to 
4,  and  2  to  3.  The  periodogram  operation  (4  to  3)  and  the  two  ACS 
estimation  operations  (1  to  2;  4  to  2)  are  irreversible  because, 
since  these  operations  are  nonlinear,  phase  information  is  lost  in 
the  generation  of  second-order  statistical  information  (power 
spectriam  and  ACS)  . 
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Several  modifications  have  been  proposed  to  improve  the  basic 
periodogram.  An  important  modification  is  the  inclusion  of  a  2-D 
window  in  the  right-hand-side  of  either  Equation  (2-99a)  or 
Equation  (2-lOla) .  In  general,  the  power  spectrum  estimate  is 
different  if  a  given  window  is  applied  to  the  ACS  estimate  or  to 
the  data  sequence.  The  numerical  differences,  however,  can  be 
small.  The  relevant  terminology  is  different  also:  a  given  window 
is  referred  to  as  a  lag  window  when  applied  to  an  ACS  estimate, 
and  as  a  data  window  when  applied  to  a  data  sequence.  Use  of  a 
window  in  either  context  leads  to  a  more  stable  power  spectrum 
estimate  (less  error  variance) ,  at  the  expense  of  less  resolution. 
When  windowed  data  is  used  to  generate  the  spectrum,  the  procedure 
is  referred  to  as  the  modified  oeriodoaram  method.  And  when  the 
windowed  ACS  is  used  to  generate  the  spectrum  (as  in  Equations  (2- 
84)  and  (2-85)  for  true  functions),  the  procedure  is  referred  to 
as  the  Blackman-Tukev  method. 

Another  approach  to  improve  the  periodogram  stability  is 
applicable  only  to  the  data-based  periodogram  definition.  Equation 
(2-lOla) .  This  approach  consists  of  three  steps:  (a)  partition 
the  given  data  sequence  into  several  shorter  segments  (with  or 
without  segment  overlap)  ,  (b)  computation  of  a  periodogram  from 
each  data  segment,  and  (c)  averaging  the  resulting  periodograms . 
A  data  window  may  be  applied  to  the  individual  sequence  segments 
(see,  for  example,  Oppenheim  and  Schafer,  1975)  .  As  with  the 
prior  modification,  enhanced  spectral  estimate  stability  is 
attained  at  the  cost  of  reduced  resolution. 


50 


3 . 0  SAMPLE  SIMULATION  RESULTS 

Several  runs  were  made  to  validate  the  analytic/software 
model,  with  emphasis  on  the  ground  clutter  component  because  it  is 
the  most  complex  segment  of  the  model.  For  that  purpose  the 
baseline  case  defined  by  Ward  (1994)  has  been  adopted,  with  the 
parameters  as  listed  in  Tables  3-1  and  3-2.  These  tables  list  the 
parameter  symbol,  description,  and  value  for  Ward's  baseline  case 
in  the  units  used  in  SSC ' s  MATLAB-based  software.  For  simplicity, 
normalized  values  are  used  for  noise,  signal,  and  jammer  powers. 

Software  validation  of  the  MATLAB-based  simulation  included 
exercised  it  with  the  values  listed  in  Tables  3-1  and  3-2  as 
inputs.  Simulation  outputs  include  the  power  ratio  parameters 
signal-to-noise  ratio  (SNR) ,  jammer- to-noise  ratio  (JNR) ,  and  CNR. 
Several  plots  are  simulation  outputs  also,  including  CNR  as  a 
function  of  azimuth,  and  true  and  estimated  2-D  spectra.  For  the 
baseline  case  in  Tables  3-1  and  3-2  the  power  ratio  parameters 
are:  0  dB  SNR,  38.1  dB  JNR,  and  46.1  dB  CNR.  Ward's  results  for 
these  same  conditions  are  (Ward,  1994) :  0  dB  SNR,  38  dB  JNR,  and 
47  dB  CNR.  The  0.9  dB  disagreement  in  the  CNR  value  is  due  to 
small  differences  between  the  two  models.  Ward  (1994)  does  not 
list  the  values  he  uses  for  the  physical  constants,  and  does  not 
provide  analytical  expressions  for  several  parameters,  such  as 
receiver  noise  power  spectral  density,  transmit  antenna  pattern, 
and  receive  column  pattern  beamforming  gain. 

Figure  3-1  presents  CNR  as  a  function  of  azimuth,  ^j.((|),0q),  for 
Ward’s  baseline  case,  with  6q  calculated  as  in  Figure  2-2.  This 

figure  corresponds  to  Figure  10  in  (Ward,  1994) ,  and  is  presented 
here  in  the  same  scale  (for  both  axes)  and  with  gridlines,  for 
ease  of  comparison.  Notice  that  the  figures  are  very  similar  in 
outline  and  scale.  The  maximum  value  of  ^j.((t),6o)  is  3  8.1  dB  for  (])  = 

0°,  whereas  the  maximum  value  extracted  visually  from  Ward's  Figure 
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10  is  approximately  38.5  dB.  This  discrepancy  is  consistent  with 
the  difference  in  total  CNR  (0.9  dB)  between  the  two  simulations. 


SYMBOL 

DESCRIPTION 

VALUE 

J 

Number  of  azimuth  axis  array  elements 
(or  sub-arrays),  which  is  the  number  of 
azimuth  channels 

18 

N 

Number  of  pulses  in  one  CPI,  which  is 
the  duration  of  the  channel  output 
vector  sequence 

18 

Je 

Number  of  elevation  axis  array 
elements  beamformed  into  one  azimuth 
channel 

4 

Pt 

Peak  transmitted  power 

200  kW 

Tu 

Uncompressed  pulse  duration 

200  |is 

Vrf 

Pulse  repetition  frequency 

300  Hz 

Radiation  frequency 

450  MHz 

Receiver  instantaneous  bandwidth 

4  MHz 

Ga(<l>.6o) 

Transmit  array  pattern 

Uniform 

Go 

Transmit  pattern  gain 

22  dB 

Ge 

Array  element  gain 

4dB 

Gb 

Array  element  backlobe  pattern  gain 

-30  dB 

Receiver  noise  power  in  each  channel 

1  (normalized) 

Fn 

Noise  figure 

3dB 

Ls 

System  losses 

4dB 

Table  3-1.  Radar  system  parameters  for  Ward  (1994)  baseline  case. 
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SYMBOL 


DESCRIPTION 


VALUE 


Hp 

Platform  altitude 

9  km 

Vp 

Platform  velocity 

50  m/s 

rc 

Range  from  platform  to  ground  clutter 
range  ring  illuminated  by  array  main 
beam 

130  km 

Y 

Aircraft  platform  crab  angle 

0  deg 

<1> 

Array  main  beam  azimuth  angle 

0  deg 

Ne 

Number  of  ground  clutter  patches  in  the 
range  ring  illuminated  by  array  main 
beam 

361 

Gs 

Number  of  moving  targets 

1 

a 

Target  model  damping  coefficient 

0.99 

Vt 

Target  radial  velocity 

33.3333  m/s 

<t>t 

Target  azimuth  angle 

0  deg 

Target  elevation  angle 

0  deg 

Target  power  in  each  channel 

1 

Gi 

Number  of  jamming  sources 

2 

‘t’i 

Jamming  sources  azimuth  angles 

25  deg;  -40  deg 

©i 

Jamming  sources  elevation  angles 

0  deg:  0  deg 

Jammer  powers  in  each  channel 

3310  dB;  3000  dB 

Table  3-2.  Scenario  parameters  for  Ward  (1994)  baseline  case. 
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CNR  as  a  function  of  azimuth  is  presented  in  Figure  3-2  for 
Ward's  baseline  case,  except  that  the  array  pattern  function, 
Ga(<t>,0o)'  ^  Dolph-Chebyshev  taper  with  30  dB  sidelobe  attenuation 

level.  For  ease  of  comparison,  the  axes  scale  is  the  same  as  in 
Figure  3-1.  As  expected,  the  CNR  reflects  the  broader  main  beam 

and  lower  sidelobes  in  relation  to  the  conditions  of  Figure  3-1. 
The  maximum  value  of  ^j,((1),0q)  is  38.1  dB  for  (|)  =  0°,  and  the  total  CNR 

is  51.6  dB.  Both  plots  have  the  same  maximum  value  because  the 
array  pattern  function  is  normalized  to  unity  at  the  main  beam 
direction  (the  constant  factor  U  in  Equation  (2-43)). 

A  plot  of  the  true  power  spectrum  for  Ward's  baseline  case  is 
presented  in  Figures  3-3  and  3-4.  Figure  3-3  is  a  3-D  plot  of  the 
logarithm  (base  10)  of  the  spectriam,  and  Figure  3-4  is  a  contour 
plot  of  the  spectrum  in  Figure  3-3.  The  spectrum  was  normalized 
so  that  its  maximum  value  is  0  dB.  This  spectrum  was  computed 
using  the  Blackman-Tukey  method,  which  consists  of  applying  the 
DFT  to  the  product  of  the  true  ACS  and  a  lag  window  (Equations  (2- 
84)  and  (2-85)).  A  Dolph-Chebyshev  window  with  60  dB  sidelobes 
was  used  to  enhance  the  spectral  features,  and  padding  with  110 
zeros  along  each  axis  was  used  to  increase  the  plot  resolution  to 
Afd  =  Afs=  1/128.  In  both  figures,  notice  the  clutter  ridge  at  a = 45° 
(corresponding  to  u  =  1),  with  mainlobe  at  1^  =  0,  fs  =  0.  Other 

relevant  parameters  calculated  by  the  program  include  main  beam 
elevation  angle,  Sq  =  -4.41°,  and  grazing  angle,  \|fg  =  3.53°.  The 
spectral  ridges  due  to  the  two  jammers,  at  fg  =  0.21  (for  <j)j^=25°)  and 
fg  = -0.32  (for  (t)|2  =  -40°),  stand  out  in  both  figures.  However,  the 
moving  target  at  =  0.334  and  fg  =  0  is  noticeable  only  in  the 

contour  plot.  Figure  2-8,  because  the  jammer  ridge  at  spatial 
frequency  fg  = -0.32  hides  it  in  Figure  3-3. 

The  plots  presented  in  Figures  3-3  and  3-4  are  based  on  the 
true  ACS  generated  as  described  in  Sections  2.1.2  through  2.1.6 
herein.  In  contrast.  Figures  3-5  and  3-6  present  3-D  and  contour 
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plots,  respectively,  of  the  logarithm  of  the  spectrum  generated 
using  the  modified  periodogram  method  applied  to  ten  independent 
zero-padded  (with  110  zeros)  pseudo-random  realizations  of  the 
channel  output  sequence  generated  as  described  in  Sections  2.1.7 
through  2.1.11.  A  Dolph-Chebyshev  data  window  with  60  dB  (power) 
sidelobes  was  applied  to  each  realization  prior  to  the  computation 
of  its  periodogram.  The  ten  periodograms  were  averaged  to  obtain 
the  periodogram  shown  in  these  two  figures.  Notice  the  similarity 
of  Figures  3-5  and  3-6  with  Figures  3-3  and  3-4,  respectively, 
which  indicates  that  the  data-based  spectrum  based  on  only  ten 
realizations  is  very  close  to  the  true  ACS-based  spectrum. 
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clutter-to-noise  ratio,  CNR  (dB)  clutter-to-noise  ratio,  CNR  (dB) 
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power  spectrum,  S(fd,fs)  (dB) 


NORMALIZED  TRUE  LOG  POWER  SPECTRUM  (Weighted  ACS) 


Figure  3-3.  Logarithm  of  the  normalized  true  power  spectrum  of 
the  array  output  for  Ward's  baseline  case  (DFT  of  the  Dolph- 

Chebyshev  weighted  ACS)  . 
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spatial  frequency,  fs 


NORMALIZED  LOG  POWER  SPECTRUM  (Modified  Periodogram) 


0^ 


spatial  frequency,  fs  0-5  0.5  doppler  frequency,  fd 


Figure  3-5.  Logarithm  of  the  normalized  power  spectrum  of  the 
array  output  for  Ward's  baseline  case  (average  of  ten  periodograms 
with  Dolph-Chehyshev  data  weighting)  . 
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spatial  frequency,  fs 


NORMALIZED  LOG  POWER  SPECTRUM  (Modified  Periodogram) 


Figure  3-6.  Contour  plot  of  the  logarithm  of  the  normalized  power 
spectrum  of  the  array  output  for  Ward's  baseline  case  (average  of 
ten  periodograms  with  Dolph-Chebyshev  data  weighting)  . 
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4 . 0  CONCLUSIONS  AND  RECOMMENDATIONS 


An  analytic  model  for  the  multichannel  signal  in  an  airborne 
surveillance  phased  array  radar  for  moving  target  detection  is 
presented  herein.  The  model  includes  target,  ground  clutter, 
broadband  interference  (jamming),  and  receiver  noise  components. 
The  ground  clutter  component  is  the  most  complex,  and  is  derived 
from  the  analytic  model  proposed  by  Jaffer  et  al .  (1991),  and 
adopted  recently  by  Ward  (1994) .  Several  diagnostic  issues  are 
discussed  also,  including  the  transformation  of  a  one-dimensional 
vector  sequence  into  a  two-dimensional  scalar  sequence,  and  the 
estimation  of  power  spectra  for  space-time  processes .  A  MATLAB- 
based  software  program  was  developed  based  on  the  analytic  model, 
and  has  been  tested  thoroughly.  The  software  model  provides  the 
means  to  generate  representative  multichannel  radar  data  in  an 
efficient  and  controlled  manner. 

The  model  discussed  herein  has  been  utilized  to  generate 
state  space  models  for  the  ground  clutter  component  in  space-time 
processes,  and  to  evaluate  the  multichannel  innovations-based 
detection  algorithm  formulated  in  Phase  I  of  the  "Multichannel 
System  Identification  and  Detection  Using  Output  Data  Techniques" 
SBIR  program.  These  results  are  presented  in  the  companion  Volume 
I  of  this  Final  Technical  Report  (Roman  and  Davis,  1996) . 

Various  enhancements  to  the  model  have  been  identified  for 
development  in  future  programs.  These  include  adding  the 
capabilities  to  model:  (a)  internal  clutter  motion;  (b)  antennae 
element  coupling;  (c)  return  from  multiple  range  rings;  (d) 
different  types  of  jammers;  and  (e)  terrain-scattered  interference 
( "hot  clutter" ) . 
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technologies  of  command,  control,  communications  and  intelligence  and  to 
transition  them  into  systems  to  meet  customer  needs.  To  achieve  this, 
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operational  capability,  readiness,  and  supportability; 
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d.  Promotes  transfer  of  technology  to  the  private  sector; 

e.  Maintains  leading  edge  technological  expertise  in  the  areas  of 
surveillance,  commimications,  command  and  control,  intelligence, 
reliability  science,  electro-magnetic  technology,  photonics,  signal 
processing,  and  computational  science. 

The  thrust  areas  of  technical  competence  include:  Surveillance, 
Communications,  Command  and  Control,  Intelligence,  Signal  Processing, 
Computer  Science  and  Technology,  Electromagnetic  Technology, 
Photonics  and  Reliability  Sciences. 


